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TEXTBOOK  MULTIGRID  EFFICIENCY  FOR  THE  INCOMPRESSIBLE 
NAVIER-STOKES  EQUATIONS:  HIGH  REYNOLDS  NUMBER  WAKES  AND 

BOUNDARY  LAYERS 

JAMES  L.  THOMAS*,  BORIS  DISKINt,  AND  ACHI  BRANDT* 

Abstract.  Textbook  multigrid  efficiencies  for  high  Reynolds  number  simulations  based  on  the  incom¬ 
pressible  Navier-Stokes  equations  are  attained  for  a  model  problem  of  flow  past  a  finite  flat  plate.  Elements 
of  the  Full  Approximation  Scheme  multigrid  algorithm,  including  distributed  relaxation,  defect  correction, 
and  boundary  treatment,  are  presented  for  the  three  main  physical  aspects  encountered:  entering  flow,  wake 
flow,  and  boundary  layer  flow.  Textbook  efficiencies,  i.e.,  reduction  of  algebraic  errors  below  discretization 
errors  in  one  full  multigrid  cycle,  are  attained  for  second  order  accurate  simulations  at  a  laminar  Reynolds 
number  of  10,000. 

Key  words,  incompressible  Navier-Stokes  equations,  textbook  multigrid  efficiency,  distributive  relax¬ 
ation,  defect-correction  iteration 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  In  the  mid-70’s,  Beam  and  Warming  [1]  presented  an  implicit  scheme  for  the  com¬ 
pressible  Navier-Stokes  equations  which  had  a  significant  impact  on  the  field  known  as  Computational  Fluid 
Dynamics  (CFD).  The  method  they  presented,  based  upon  a  spatial  factoring  of  the  implicit  equations  in 
delta  form,  used  alternating  tridiagonal  line  relaxations  to  solve  high  Reynolds  number  viscous  simulations. 
This  method  proved  to  be  much  more  efficient  than  other  approaches.  The  basic  methodology  is  still  widely 
used  and  has  been  extended  to  very  general  applications  across  the  Mach  number  range,  forming  the  foun¬ 
dation  for  many  general  purpose  solvers  worldwide,  among  them  ARC3D  [8]  and  CFL3D  [7]  at  the  NASA 
Ames  and  Langley  Research  Centers,  respectively.  This  seminal  contribution  of  Beam  and  Warming  was  a 
critical  building  block  to  the  acceptance  of  CFD  using  Reynolds- Averaged  Navier-Stokes  (RANS)  solvers  by 
the  aircraft  industry.  Today,  computational  methods  for  the  cruise  shapes  of  transport  aircraft,  designed 
to  minimize  viscous  and  shock  wave  losses  at  transonic  speeds,  are  reasonably  well  in  hand.  Simulations  of 
off-design  performance,  involving  unsteady  separated  and  vortical  flows  with  stronger  shock  waves,  require 
significantly  greater  computing  resources;  this  requirement  limits  further  inroads  into  the  design  process  with 
CFD. 

As  a  typical  example  of  current  RANS  capability,  the  CFL3D  code  is  based  on  the  spatially-factored 
scheme  of  Beam  and  Warming  and  uses  multigrid  to  accelerate  convergence  to  steady  state;  using  alternating- 
line  implicit  block  5x5  matrix  solutions,  approximately  200  updates  are  required  to  converge  the  lift  and 
drag  to  one  percent  of  their  final  values  for  wing-body  geometries  near  transonic  cruise  conditions.  Complex 
geometry  and  complex  physics  simulations  generally  require  many  more  residual  evaluations  to  converge,  and 
sometimes  convergence  cannot  be  attained.  Now,  it  is  well-known  for  fully  elliptic  problems  that  solutions 
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can  be  attained  using  a  full  multigrid  (FMG)  process  in  far  fewer,  on  the  order  of  2-4,  residual  evaluations. 
Optimal  convergence  is  defined  by  Brandt  [2,  3,  4]  as  textbook  multigrid  efficiency  (TME),  meaning  the 
solutions  to  the  governing  system  of  equations  are  attained  in  a  computational  work  which  is  a  small  (less 
than  10)  multiple  of  the  operation  count  in  the  discretized  system  of  equations.  Thus,  there  is  a  potential 
gain  of  more  than  two  orders  of  magnitude  in  operation  count  reduction  if  TME  could  be  attained  for  the 
RANS  equation  sets.  The  principal  difficulty  stems  from  the  fact  that  the  RANS  equation  sets  are  a  system 
of  coupled  nonlinear  equations  which  are  not,  even  for  subsonic  Mach  numbers,  fully  elliptic,  but  contain 
hyperbolic  factors.  Brandt  [4]  has  summarized  the  progress  and  remaining  barriers  to  achieving  TME  for 
the  equations  of  fluid  dynamics. 

The  purpose  of  this  paper  is  to  present  a  multigrid  method  which  attains  textbook  efficiencies  for  one 
of  the  most  basic  simulations  encountered  in  fluid  dynamics  -  the  incompressible  viscous  flow  past  a  finite 
flat  plate  at  high  Reynolds  number.  The  flow,  although  relatively  simple,  contains  several  basic  elements  of 
the  barriers  to  be  overcome  in  extending  textbook  efficiencies  to  the  compressible  RANS  equations,  namely 
entering  flows,  far  wake  flows,  and  boundary  layers.  A  central  element  of  the  multigrid  method  presented 
is  the  decomposition  through  distributed  relaxation  [3]  of  the  the  system  of  equations  into  separate,  usually 
scalar,  factors  that  can  be  treated  optimally,  i.e.,  through  marching  for  the  hyperbolic  factors  and  through 
multigrid  for  the  elliptic  factors.  Although  we  restrict  ourselves  to  incompressible  flow,  the  procedures  carry 
over  directly  to  the  compressible  flow  case,  at  least  for  subcritical  flow  [3,  4,  10]. 

2.  Governing  Equations.  The  equations  considered  here  are  the  steady,  incompressible  Navier-Stokes 
equations  in  nonconservative  form,  i.e.,  two  momentum  equations  and  the  continuity  equation, 


r(q)  s  Cq  =  0, 


expressed  in  terms  of  primitive  (velocities  and  pressure)  variables  q  =  (u,v,p)T,  where 


(2.1) 


r  Qu  o  ex  i 

£  =  I  0  Qu  dy  I  .  (2.2) 

'  8X  dy  O' 

The  operator  Qu  represents  convection  and  diffusion  effects  as 


Qu  =  Q~  i/A,  (2.3) 

where  Q  =  udx  +  vdyi  the  Laplacian  operator  is  A  =  dxx  +  dyy ,  and  the  kinematic  viscosity  is  v  =  1  /Re, 
where  Re  is  Reynolds  number.  Extensions  to  conservation  law  form  for  the  momentum  equations  and  to 
inclusion  of  the  energy  equation  are  possible,  but  not  considered  here. 

The  determinant  of  the  matrix  of  operators, 


\C\  =  -Qt/A,  (2.4) 

corresponds  to  an  elliptic  factor,  represented  by  the  Laplacian,  and  a  convection-diffusion  factor,  gener¬ 
ally  recognized  as  the  convection  and  diffusion  of  vorticity  along  a  streamline.  For  high  Reynolds  number 
simulations,  there  are  two  important  scales:  the  viscous  scales  in  the  thin  viscous  layers  near  bodies  and 
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in  their  wakes  and  the  inviscid  scales,  which  predominate  over  most  of  the  flow  field.  For  the  numerical 
calculations  below,  the  thin-layer  approximation,  in  which  only  the  viscous  terms  associated  with  variations 
in  the  coordinate  normal  to  the  body  are  retained,  is  used. 

3.  Multigrid  Method.  The  present  approach  uses  a  full  multigrid  (FMG)  algorithm  [2,  3],  proceeding 
from  the  coarsest  grid  to  finer  grids.  The  solution  is  interpolated  from  the  current  grid  to  the  next  finer 
grid.  The  goal  of  the  algorithm  is  fast  reduction  of  the  algebraic  errors  below  the  discretization  errors  on  a 
given  grid,  before  moving  to  the  next  finer  grid.  The  algebraic  errors  of  the  discrete  equations  on  a  given 
grid  are  reduced  through  a  Full  Approximation  Scheme  (FAS)  [3]  multigrid  scheme,  in  which  corrections  to 
the  nonlinear  equations  are  obtained  from  coarser  grid  solutions.  The  scheme  is  described  below  by  means 
of  a  two-grid  notation,  in  which  the  fine  grid  is  denoted  by  superscript  h  and  the  coarse  grid  by  superscript 
2  h. 

The  steady-state  residual  operator  to  be  solved  on  the  fine  grid  is  the  discrete  version  of  Eq.  (2.1), 


r^(q^)  =  0. 

The  initial  fine-grid  approximation  qh  is  prolonged  from  the  coarse-grid  solution  q2/l,  as 


(3.1) 


qh  <r-  Vq2h.  (3.2) 

where  V  denotes  a  prolongation  operator.  After  relaxation(s)  of  the  fine-grid  operator  to  obtain  an  approx¬ 
imation  qh,  the  coarse-grid  equation  at  level  2 h  to  be  solved  for  a  correction  to  the  fine  grid  is 


r2h(q2h)  =  v2h{Uqh)  -  ftrh(qh),  (3.3) 

where  1Z  denotes  a  restriction  operator  for  transfer  of  information  to  the  coarser  grid  and  the  tilde  superscript 
denotes  a  most  recently  available  value.  This  coarse  grid  equation  is  then  solved  by  some  iterative  method 
(or  directly  if  the  grid  is  coarse  enough).  The  correction  from  the  coarse  grid  (grid  2 h)  is  prolonged  to  the 
finer  grid  as 


qh  <-  qh  +  V(q2h  -  7Zqh).  (3.4) 

The  restrictions  1Z  used  here  are  volume-weighted  for  the  continuity  equations;  for  the  momentum  equations, 
the  coarser  cell  values  are  found  by  volume- weighted  restrictions  in  the  direction  parallel  to  the  cell  interface 
along  with  full- weighted  restrictions  in  the  orthogonal  direction,  i.e.,  for  the  y-~ momentum  equation,  volume- 
weighted  horizontally  and  full-weighted  vertically.  The  prolongations  V  are  bicubic  interpolations  from 
coarser  meshes  for  both  the  solution  and  the  correction  ,  although  results  with  linear  interpolations  were 
nearly  identical.  The  FAS  cycle  described  above  is  used  extensively  in  current  Euler  and  Navier-Stokes 
solvers.  The  algorithm  is  critically  dependent  on  the  choice  of  relaxation  operator;  distributed  relaxation  is 
used  here  as  described  subsequently. 

The  coarse-grid  equations  are  themselves  solved  with  7  cycles  of  the  algorithm  applied  recursively,  where 
7=1  would  correspond  to  a  V-cycle  and  7  =  2  to  a  W-cycle;  the  number  of  relaxations  on  the  downward 
and  upward  legs  of  the  cycle  are  denoted  as  (^1,^2).  We  use  here  (i/i,z/2)  —  (2,1)  and  a  variant  of  the 
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FV(v1tv2)  Cycle 


FlG.  3.1.  Schematic  of  the  FV-cycle  for  4-level  multigrid  where  vq  denotes  the  number  of  relaxations  on  the  coarsest  mesh 

(Cl8h  ). 

V-cycle,  termed  an  FV-cycle,  in  which  the  initial  approximation  to  the  correction  on  the  2h  grid  is  obtained 
through  a  FMG  process.  The  cycle  is  sketched  in  Fig.  3.1;  the  amount  of  additional  computational  work 
compared  to  a  standard  V-cycle  is  small,  in  the  ratio  of  8/7  in  the  limit  of  an  infinite  number  of  levels  in  two 
dimensions.  For  the  simulations  here,  six  levels  were  used  wherever  possible.  The  notation  FMG-n  denotes 
an  FMG  cycle  with  n  FV(2,1)  cycles  at  each  level. 

4.  Distributed  Relaxation.  Away  from  boundaries,  the  correction  Sq  to  the  current  approximation 
q ,  introduced  at  the  stage  of  distributed  relaxation,  [2,  3]  is  calculated  from 

L<5q  =  -r(q),  (4.1) 

where  L  is  a  principal  linearization  of  £,  in  which  the  coefficients  u  and  v  in  Eq.  (2.3)  are  evaluated  from 
the  current  approximation  and  fixed  throughout  the  relaxation.  Note  this  is  not  a  Newton  linearization; 
only  the  principal  terms  at  the  viscous  and  inviscid  scales  are  retained.  The  distributed  relaxation  method 
replaces  6q  by  M5w  so  that  the  resulting  matrix  LM  becomes  a  diagonal  or  lower  triangular  matrix,  as 

LMiSw  =-r(q).  (4.2) 

The  diagonal  elements  of  LM  are  composed  ideally  of  the  separate  factors  of  the  determinant  of  the  matrix 
L  and  represent  the  elliptic  or  hyperbolic  features  of  the  equation.  For  incompressible  flow,  the  distribution 
matrix  M  can  take  on  a  particularly  simple  form,  as  determined  by  the  cofactors  of  the  third  row  of  L 
divided  by  their  common  factor,  as 

'  1  0  -Qxm 

M  =  0  1  -By  ,  (4.3) 

_  0  0  Qv 
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yielding 


r  Qv  o  o  i 

LM  =  I  0  Qu  0  1.  (4.4) 

dx  dx  -A 

The  determinant  of  the  operator  matrix  LM, 

|LM|  =  -e#A,  (4.5) 

corresponds  to  an  elliptic  factor  and  two  convection-diffusion  factors;  the  additional  term  over  Eq.  (2.4),  Qu, 
indicates  that  as  a  set  of  new  variables,  Sw  would  generally  need  additional  boundary  conditions  all  around 
the  boundary  (or,  just  at  inflow  in  the  case  v  =  0).  Brandt  termed  the  variables  Sw  as  “ghost  variables,” 
since  they  need  not  explicitly  appear  in  the  calculations;  here,  they  do  appear  in  the  calculations,  although, 
as  with  the  original  intent,  the  boundary  conditions  are  derived  from  the  original  primitive  variables.  The 
equations  to  solve  for  the  ghost  variables  are  given  explicitly  as 


QiySwi  -  -n, 

Qv6w2  =  -r2, 

ASws  =  +r3  +  dxSwi  +  dySw2- 

(4.6) 

Near  boundaries,  the  general  approach,  [3,  4]  would  be  to  relax  the  governing  equations  directly,  since 
the  equations  do  not  necessarily  decouple  near  boundaries  as  they  do  in  the  interior  of  the  domain.  One 
can  make  more  general,  but  possibly  slowly  converging,  relaxations,  such  as  Kaczmarcz  relaxation,  in  this 
region.  This  will  not  affect  the  overall  complexity,  because  the  number  of  boundary  points  is  negligible  in 
comparison  to  the  number  of  interior  points.  Here,  however,  we  use  an  approach  which  applies  the  interior 
distributive  relaxation  operator  also  at  the  boundaries,  inferring  boundary  conditions  for  the  ghost  variables 
based  on  the  boundary  conditions  of  the  governing  equations.  The  cost  is  that  the  correction  equations, 
Eq.  (4.6),  no  longer  assume  a  triangular  form,  requiring  a  block  matrix  solution  at  the  boundaries  rather 
than  the  scalar  solutions  attained  away  from  the  boundary.  Assuming  linearized  flow,  the  appropriate  ghost 
variable  boundary  conditions  at  the  differential  level  are  derived  for  inviscid  inflow  and  outflow  in  Appendix 
I  and  tangency  in  Appendix  II  .  These  boundary  conditions  are  implemented  discretely  at  the  corresponding 
boundaries.  At  the  no-slip  boundary,  the  corresponding  discrete  boundary  conditions  for  the  ghost  variables 
are  constructed  in  Appendix  III.  The  procedure  is  effective  for  the  simulations  considered  here;  details  are 
given  in  subsequent  sections. 

5.  Defect  Correction  Relaxation.  Since  Eq.  (4.2)  is  written  in  delta  form,  it  is  natural  to  consider 
defect  correction  for  the  update,  namely  a  lower-order  discretization  of  the  left  side  of  Eq.  (4.6)  in  order 
to  simplify  the  construction  and  reduce  the  bandwidth  of  the  implicit  operator.  Here,  we  use  a  first-order 
upwind  discretization  for  the  convective  part  of  the  convection-diffusion  operator,  Qv,  in  Eq.  (4.6).  The 
distributed  relaxation  operator  can  thus  be  written  as 


[L  M]d  Sw  =  -rt, 


(5.1) 
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Fig.  6.1.  Variable  description  for  a  grid  of  Jx K  —  NxxNy  points. 


where  the  subscripts  t  and  d  denote  some  desired  “target”  and  “driver”  schemes  on  the  right  and  left  sides, 
respectively,  of  the  equation. 

For  hyperbolic  equations,  the  initial  convergence  of  defect  correction  may  be  slow  for  certain,  not  nec¬ 
essarily  high,  frequencies  [9,  10,  6].  For  a  target  second-order  upwind-biased  discretization  corresponding 
to  k  =  0,  defined  subsequently,  the  asymptotic  convergence  rate  is  approximately  0.5  per  defect-correction 
iteration.  Thus,  it  is  well-matched  with  the  convergence  rate  of  0.5  per  relaxation  expected  for  the  elliptic 
parts  of  the  operator  with  Gauss-Seidel  relaxation. 

Defect  correction  is  implemented  in  the  multigrid  algorithm  as  follows:  any  discrete  evaluations  of 
the  residuals  of  Eq.  (3.1)  (including  residuals  transferred  to  the  coarse  mesh)  are  done  with  the  target 
discretization  and  any  updates  via  distributed  relaxation  are  done  with  the  driver  operator,  which  is  first- 
order  upwind  for  convection.  This  is  similar  to  the  “double-discretization”  approach  of  Brandt  [3]  in  practice, 
except  that  the  target  residual  is  evaluated  on  all  of  the  meshes,  including  the  finest  mesh. 

6.  Numerical  Discretization. 

6.1.  Spatial  Discretization.  The  st agger ed-gr id  discretization  used  here,  as  shown  in  Fig.  6.1,  is  usual: 
p  defined  at  the  cell-centers  of  the  grid,  u  defined  at  the  cell  interfaces  tangent  to  the  y—  or  A:— direction,  and 
v  defined  at  the  cell  interfaces  tangent  to  the  a;—  or  j -direction.  Additional  values  of  v  and  p  are  defined 
at  inflow  and  outflow  boundaries  in  order  to  accommodate  boundary  conditions,  defined  subsequently.  The 
discrete  scheme  with  such  a  staggered-grid  arrangement  of  variables  can  be  described  as 

'  Qt  o  an 

LVs  0  Q’i  a £  qh  =  0,  (6.1) 

.  a£  a£  o  . 
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where  Qj  and  dy  are  generally  distance-h  central  differences  on  the  staggered  grid.  The  operator  is 
composed  of  convection  and  diffusion  elements,  analagous  to  Eq.  (2.3);  the  diffusion  elements  are  treated 
with  central  differencing, 


(Q  u)hu  =  -  --p’fc+1  “ 

(  yy  hk  ( hy)jtk[  (hy)jM1/2 


^j,k  'U'jyk  —  lj 

(hy)j,k- 1/2 


(6.2) 


where  hy  denotes  grid  spacing  in  the  y  direction. 

The  discrete  convection  operator  Qh  is  upwind-biased,  of  either  the  standard  upwind  differencing  (SUD) 
type  or  the  narrow  upwind  differencing  (NUD)  type.  The  operator  can  be  defined  on  a  uniform  grid  in  terms 
of  translation  operators  Tfm  and  Tfc±m,  (T^mujik  =  uj±m,k).  The  SUD  scheme  can  be  defined  as 


Qh  =  p.D(T;an{u))  +  M  D(T^v)) 

flx  hy 

where  hx  is  the  grid  spacing  in  the  x  direction,  the  sign  function  sgn  is  defined  as 

(  +1  ifx>0, 
sgn(x)  —  <  -1  ifz<0, 

[  0  otherwise, 

and  D  is  defined  as 


(6.3) 


D(z)  =  C-2Z  H-  C—\Z  1+Co  +  Ci2?+1. 
The  NUD  scheme  can  be  defined  as  below  for  j^j-  >  , 


Qh  =  ( N  _  M )2P(T'?5n(u))  -f  —  D(T^gn^ xS9n^) 


M; 


' h . 


and  as  below  for  <C 


(6.4) 


Qh 


^ D(T • 
hx 


,sgn{u)rj.sgn(v _|_  (.M  _  1h! 

hy  hx 


(6.5) 


For  uniform  meshes  or  meshes  in  which  the  stretching  ratio  is  ft  =  1  4-  0(h),  k— schemes  of  at  least  second 
order  accuracy  (SUD-2  and  NUD-2)  are  defined  for  «  €  [— 1, 1]  as 


{c_2,c_i,co,ci}  =  2+1^{l  -  «i  3/6-  5,  3(1  -  k),  1  +  k} 

and  third-order  accuracy  (SUD-3)  is  attained  for  k  =  1/3  with  uniform  meshes  or  meshes  in  which  the 
stretching  ratio  is  ft  —  1  -F  0(h2).  On  stretched  grids,  the  reference  meshsize,  hi- 1/2,  appearing  in  (the 
denominator  of)  the  discrete  one-dimensional  convection  operator,  D(Ti)/hi- 1/2,  is  a  meshsize  upstream  of 
the  i-th  node  where  the  discrete  operator  is  defined.  The  coefficients  for  the  first-order  upwind  schemes 
(SUD-1  and  NUD-1)  are 


{c_2,c-i,co,Ci}  =  {0,  -1,  1,  0}. 
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Table  7.1 

Errors  in  u  with  the  FMG-1  cycle  for  entering  flow  using  a  second  order  accurate  discretization  of  the  continuity  equation; 


Scheme 

h 

IMI  :  u 

IMI/IMI  :u 

SUD-l 

1/16 

0.115556x10° 

0.019 

SUD-1 

1/32 

0.664116xl0_1 

0.008 

SUD-l 

1/64 

0.357011xl0_1 

0.006 

SUD-1 

1/128 

0.185119xl0_1 

0.002 

NUD-1 

1/16 

0.476075xl0-1 

0.007 

NUD-1 

1/32 

0.246260x10“ 1 

0.008 

NUD-1 

1/64 

0.125445xl0_1 

0.006 

NUD-1 

1/128 

0.633386xl0-2 

0.003 

SUD-2 

1/16 

0. 68900 lxlO-2 

0.024 

SUD-2 

1/32 

0.154126xl0-2 

0.039 

SUD-2 

1/64 

0.368421xl0-3 

0.034 

SUD-2 

1/128 

0.905679xl0-4 

0.026 

NUD-2 

1/16 

0.251242xl0-2 

0.128 

NUD-2 

1/32 

0.637956xl0-3 

0.046 

NUD-2 

1/64 

0.159458xl0-3 

0.046 

NUD-2 

1/128 

0.397594xl0-4 

0.047 

Table  7.2 

Errors  in  u  and  v  for  entering  flow  with  a  fourth  order  accurate  discretization  of  the  continuity  equation;  t=0.5. 


Scheme 

h 

IMI 

IMh« 

SUD-3 

1/32 

0.325327xl0-3 

0.182866xl0-3 

SUD-3 

1/64 

0.425745xl0“4 

0.231228xl0~4 

SUD-3 

1/128 

0.547187xl0-5 

0.292787xl0-5 

NUD-2 

1/32 

0.121481xl0-3 

0.591146xl0-4 

NUD-2 

1/64 

0.151229xl0-4 

0.724727xl0-5 

NUD-2 

1/128 

0.186856xl0-5 

0.885510x10-° 

6.2.  Gauss-Seidel  Line  Relaxation.  The  equations  for  <5w  are  relaxed  with  a  line-y  Gauss-Seidel 
algorithm  marching  from  the  inflow  to  the  outflow  boundary.  The  correction  equations  for  Sw  are  solved 
implicitly  because  of  the  highly  stretched  mesh  used  for  the  viscous  calculations.  Since  the  thin-layer  ap¬ 
proximation  is  made  for  the  viscous  terms,  the  convective  operator  is  first-order  upwind,  and  there  is  no 
streamwise  reversed  flow,  the  Sw±  and  Sw2  correction  (driver)  equations  of  Eq.  (4.6),  corresponding  to  the 
linearized  momentum  equations  at  given  pressure,  are  solved  exactly.  The  line-y  solutions  require  only  in¬ 
versions  of  tridiagonal  (rather  than  block-tridiagonal)  matrices,  since  the  equations  for  Sw  form  a  lower 
triangular  set  except  near  the  boundaries.  The  treatment  at  boundaries  requires  special  consideration  as 
discussed  subsequently.  Note  that  for  the  NUD  schemes  in  inviscid  flow  with  the  tridiagonal 

equations  for  Sw±  and  Sw2  reduce  to  diagonal  equations. 
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•© -  Algebraic/Discretization 

-b -  Residual 


-e -  Algebraic/Discretization 

■e -  Residual 


(b)  Third  order  accurate  solution. 

Fig.  7.1.  Maximum  residual  and  algebraic- to- discretization  errors  in  u  versus  multigrid  cycle  for  the  three  finest  grids 
with  the  NUD-2  scheme. 


7.  Entering  Flow  Simulation.  The  flow  field  upstream  of  an  external  aerodynamic  simulation  is 
basically  inviscid.  Brandt  and  Yavneh[5]  considered  multigrid  solutions  of  such  flows  and  showed  the  accuracy 
of  the  coarse  grid  correction  to  be  critically  dependent  on  the  alignment  of  the  flow  relative  to  the  mesh.  Their 
numerical  results  indicated  the  necessity  of  W-cycles  to  converge  the  algebraic  errors  below  discretization 
errors  in  the  FMG-1  cycle.  We  revisit  these  simulations  below  with  slightly  different  boundary  conditions 
and  show  that  the  FMG-1  cycle  with  the  use  of  FV-cycles  is  sufficient.  The  computations  were  done  for 
a  square  domain  with  periodicity  in  the  ^/-direction  on  a  uniform  mesh.  Inflow  boundary  conditions  were 
specified  velocities  as 
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(7.1) 


u{0,y )  =  1  +  0.5cos(27ry), 

w(0,J/)  =  M°>J/)> 

with  constant  pressure  at  the  outflow  boundary.  The  tangent  of  the  angle  of  the  flow  relative  to  the  grid 
is  t  =  0.5,  corresponding  to  the  maximum  value  studied  by  Brandt  and  Yavneh[5].  The  exact  solution 
corresponds  to  convection  along  a  streamline  at  constant  pressure, 

u(x,  y)  =  1  +  0.5cos(27r(y  -  tx)), 

v(x,y)  =  tu(x,y).  (7-2) 

The  boundary  conditions  for  the  correction  equations  are  implemented  by  applying  the  distributed 
relaxation  equations  <5q  =  M<5w  at  the  boundary  along  with  a  Dirichlet  condition  for  {Sw)3  at  inflow.  The 
resulting  discrete  boundary  conditions  at  x  —  0  are 


5w2  =  0, 

Sws  =  0, 

Swi  =  <9£(<5w3).  (7-3) 

This  boundary  condition  is  the  discrete  equivalent  to  the  original  problem  statement  for  the  constant  coef¬ 
ficient  problem.  This  boundary  condition  couples  the  Sw i  and  6W3  equations  together  at  the  line  of  cells 
adjacent  to  the  inflow  boundary,  necessitating  a  block  2x2  block  matrix  solution  procedure;  away  from  this 
first  line,  the  equations  retain  the  triangular  form  of  Eq.  (4.6)  and  can  be  solved  as  scalar  equations.  The 
downstream  boundary  condition  is  implemented  by  solving  for  S1U3  at  the  last  interior  column  of  cells  simul¬ 
taneously  with  Sw 3  at  the  outflow  column,  again  necessitating  a  2x2  block  matrix  tridiagonal  solution.  After 
sweeping  through  the  domain,  all  of  the  momentum  equation  residuals  are  zero  in  the  constant  coefficient 
case;  this  local  block  matrix  coupling  at  either  boundary  eliminates  the  need  for  the  extra  sweep  of  the 
residual  equation  advocated  by  Brandt  and  Yavneh[5].  The  residuals  remain  non-zero  in  the  general  case 
because  of  subprincipal  terms  and  are  restricted  to  the  coarse  grids.  Enforcing  periodicity  in  the  y-direction 
in  the  tridiagonal  solver  eliminates  the  need  to  consider  any  special  boundary  conditions  in  that  direction. 
Special  forms  for  the  spatial  discretization  of  the  convective  operator  in  Eq.  (6.1)  at  inflow  and  outflow  are 
given  in  Appendix  TV. 

The  L2— norms  of  the  discretization  errors  in  u  after  complete  convergence  and  the  ratios  of  the  L2-norm 
of  the  algebraic  errors  divided  by  the  L2-norm  of  the  discretization  errors  after  one  cycle  are  shown  in 
Table  7.1  for  various  grid  sizes  and  orders  of  accuracy.  The  algebraic  errors  are  reduced  substantially  below 
the  discretization  errors  in  one  cycle.  The  error  norms  indicate  a  first  order  accuracy  for  SUD-1  and  NUD-1, 
and  second  order  for  SUD-2  and  NUD-2,  as  expected. 

At  this  flow  angle,  t  =  0.5,  the  NUD-2  scheme  exhibits  third-order  accuracy  for  the  linearized  convection 
problem  but  does  not  for  the  full  Euler  equations  because  second  order  accurate  discretizations  are  used  for 
the  continuity  equation,  for  the  pressure  terms  in  the  momentum  equation,  and  for  the  reconstruction  of 
the  flow  at  an  interface.  To  remedy  this,  these  discretizations  were  improved  to  fourth  order  accuracy;  the 
corresponding  results  shown  in  Table  7.2  for  both  the  SUD-3  and  NUD-2  schemes  now  exhibit  third-order 
accuracy  in  u  and  v. 
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(x,y)=(0,l) 


Fig.  8.1.  Grid  used  for  the  wake  and  finite  flat  plate  simulation. 

Table  8.1 

Computed  values  of  centerline  velocity  at  x  =  1.5  for  the  wake  simulation;  SUD-2  scheme ;  «  =  0;  Wd  —  0.5;  Re  —  10,000. 


iV^XNy 

u 

(FMG-10) 

u 

(FMG-1) 

IMI/IMI 

(FMG-1)  . 

49x25 

0.730529 

0.730585 

0.00445 

97  x  49 

0.740382 

0.740412 

0.01135 

193  x  97 

0.742367 

0.742385 

0.02672 

The  reduction  of  the  maximum  residual  and  the  algebraic-to-discretization  errors  over  4  cycles  for  the 
three  finest  grids  in  the  calculation  are  shown  in  Fig.  7.1  for  the  NUD-2  scheme  with  second  and  fourth  order 
accurate  discretizations  of  the  continuity  equation.  For  second  order  accuracy,  the  residual  and  algebraic- 
to-discretization  errors  are  reduced  four  orders  of  magnitude  over  the  4  cycles,  close  to  the  theoretical  limit 
expected  for  elliptic  equations  of  (0.5)3  =  0.125  reduction  per  FV(2,1)  cycle.  The  convergence  for  the  third 
order  accurate  results  deteriorate  somewhat  to  three  orders  of  magnitude  over  the  four  cycles  but  is  still 
quite  reasonable  considering  that  defect  correction  with  a  first-order  driver  operator  is  being  used.  Further 
improvements  could  be  made  by  additional  sweeps  or  by  a  predictor-corrector  sequence  of  the  momentum 
equations  only,  since  the  deficiency  resides  with  the  first-order  accuracy  in  the  driver  operator  for  convection. 

8.  Wake  Flow  Simulation.  The  wake  and  the  finite  flat  plate  simulation  to  follow  were  computed  for 
the  computational  domain  shown  in  Fig.  8.1  at  a  Re=10,000  based  on  the  height  of  the  channel.  The  grid 
was  stretched  in  the  y -direction  with  a  stretching  factor  on  a  specified  mesh  defined  as 

A  =  (hy)j,k+l/{hy)j,k 

corresponding  to  ( Ny)0  grid  points  in  the  vertical  direction.  The  stretching  ratio  on  all  other  meshes  is 

Freestream  pressure  is  specified  at  the  outflow  boundary;  a  wake  deficit  was  prescribed  at  the  inflow  boundary, 
x  —  0,  according  to 
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Cycles 


(a)  1/2 —norm  of  the  residual. 


Cycles 


(b)  Algebraic-to-discretization  errors  in  mass  flow. 

Fig.  8.2.  Wake  simulation  convergence  using  the  FMG-5  cycle;  Wd  —  0.5/  SUD-2  scheme;  k  —  0. 


— Re  v2 

u(0,y)  =  1  -u^exp( — v(0,y)  =  0, 

where  Wd  =  0.5.  The  mass  flow  is  defined  as  the  integral  of  velocity  at  constant  x\  the  exact  value  is 
0.9911377307.  The  boundary  condition  treatments  at  inflow  and  outflow  are  the  same  as  those  for  the 


Fig.  9.1.  Convergence  of  the  Cf  values  at  x  —  1.5  with  nominal  grid  spacing  for  two  stretching  ratios ;  SUD-2  scheme; 
■  0;  Re  —  10, 000. 


entering  flow  discussed  previously.  Symmetry  conditions  are  applied  at  y  =  0  for  both  q  and  Sw.  A 
tangency  condition,  v  =  0,  was  applied  at  y  =  1;  applying  the  distribution  operator  at  this  point  with 
simple  reflection  for  u  across  the  boundary  indicates  that  a  Neumann  condition  can  be  applied  to  6ws  at 
the  boundary,  as  shown  in  Appendix  II,  along  with  reflection  for  8w\  and  a  Dirichlet  condition  for  5wz ,  if 
needed. 

The  finest  grid  considered  was  A^xNy  —  193x97  with  /3q  =  1.03  corresponding  to  ( Ny)o  =  97.  In 
addition  to  residuals,  the  centerline  velocity  (obtained  by  second-order  extrapolation)  and  the  mass  flow 
were  monitored  at  x  =  1.5,  a  location  midway  in  the  domain,  as  a  measure  of  spatial  convergence.  The 
reductions  of  the  maximum  residual  and  the  algebraic-discretization  errors  in  mass  flow  for  all  the  grids  in 
an  FMG-5  process  are  shown  in  Fig.  8.2  using  the  SUD-2  scheme.  For  each  of  the  meshes,  the  residual  is 
reduced  4-5  orders  of  magnitude  and  the  algebraic  errors  are  reduced  far  below  discretization  errors.  The 
centerline  velocities  for  the  three  finest  grids,  Table  8.1,  demonstrate  second  order  accuracy  with  algebraic 
errors  reduced  below  discretization  errors  using  the  FMG-1  cycle.  The  reference  centerline  velocity  was 
obtained  by  second  order  Richardson  extrapolation. 

Although  not  shown,  parameter  variations  in  Wd  were  made  which  indicated  the  results  were  not  sensitive 
to  Wd  over  the  range  investigated,  0  to  0.9.  This  is  in  contrast  to  an  earlier  application, [10]  in  which  the 
ghost  variable  equations  were  solved  with  a  correction  scheme  (CS)  multigrid.  Those  results  deteriorated  for 
high  values  of  Wd ,  emphasizing  the  advantage  of  applying  the  FAS  multigrid  scheme  to  the  whole  nonlinear 
system  of  equations.  For  linear  equations,  the  performance  of  FAS  multigrid  is  the  same  as  CS  multigrid  . 

9.  Flat  Plate  Boundary  Layer  Simulation.  For  the  flat  plate  simulation,  no-slip  conditions  are 
prescribed  from  x  =  1  to  x  =  2  along  the  lower  boundary  and  symmetry  conditions  upstream  and  downstream 
of  those  points;  a  wake  profile  develops  downstream  of  the  trailing  edges,  x  =  2.  The  inflow  and  outflow 
conditions  are  prescribed  freestream  velocities  (u^  =  l,Uoo  =  0)  and  pressures,  respectively.  The  discrete 
velocities  adjacent  to  the  plate  for  y  <  0  are  required  to  satisfy  the  no-slip  condition  at  the  plate,  i.e. 
u(x,-hy/ 2)  =  ~u{x,hy/2)-,v(x,  -hy)  =  -v{x,hy).  The  distributive  relaxation  equations  applied  at  the 
boundary  are  shown  in  Appendix  III. 
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Table  9.1 

Computed  values  of  total  drag  for  the  finite  flat  plate  simulation j  SUD-2  scheme /  k  —  0/  Re  10, 000,  (3q  1.03. 


NXX  Ny 

CD 

(FMG-10) 

Cd 

(FMG-1) 

IMI/IMI :  CD 
(FMG-1) 

49x25 

0.011552 

0.011753 

0.0784 

97x49 

0.013492 

0.013412 

0.1284 

193  x  97 

0.013961 

0.014051 

0.5760 

Fig.  9.2.  Pressures  (y  =  0)  for  the  finite  flate  plate ;  Re  —  10, 000. 


The  spatial  convergence  of  the  local  skin  friction  Cf  midway  down  the  plate  versus  the  nominal  grid 
spacing  for  two  families  of  meshes  for  two  stretching  ratios  is  shown  in  Fig.  9.1,  where 


Cf  =  2v{dyu)/u^x>. 

The  two  finest  grids  in  each  family  are  289x145  and  193x97.  Second  order  accuracy  is  evident;  the  results 
with  higher  stretching  ratio  are  slightly  more  accurate  on  coarser  grids.  The  results  converge  to  a  value 
approximately  five  percent  higher  than  the  Blasius  value,  Cf  =  0.664/ y/Rex  =  0.00939,  where  x  denotes 
distance  from  the  leading  edge,  because  of  the  presence  of  a  favorable  pressure  gradient  (accelerating  flow) 
over  most  of  the  plate,  as  shown  in  Fig.  9.2.  Convergence  of  the  L2-norm  of  the  residual  and  estimated 
algebraic-to- discretization  errors  in  total  drag  Cd  are  shown  in  Fig.  9.3.  The  total  drag  is  defined  as 

Cd  —  2Cf(x*  —  l)+f  Cfdx , 

J x* 

where  the  Cf  behavior  ahead  of  x*  —  1.25  is  assumed  to  be  an  inverse  square  root  behavior  in  distance  from 
the  leading  edge,  as  occurs  with  the  Blasius  solution.  The  infinite-grid  result  is  extrapolated  using  the  two 
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(a)  L2-norm  of  the  residual. 


(b)  Algebraic-to-discretization  errors  in  Cd- 
Fig.  9.3.  Errors  per  cycle  using  the  FMG-5  cycle;  fio  =  1.03;  SUD-2  schemes  =  0;  Re  =  10,000. 


finest  grids.  Both  the  residual  and  algebraic-to-discretization  errors  are  reduced  nearly  four  orders  of  mag¬ 
nitude  over  five  cycles  for  the  four  finest  grids,  close  to  the  convergence  expected  for  elliptic  equations.  The 
Cd  values  on  the  three  finest  meshes  are  given  in  Table  9.1,  confirming  that  the  algebraic-to-discretization 
errors  are  reduced  below  unity  in  a  single  cycle.  The  values  extrapolate  to  a  slightly  larger  value  than  the 
Blasius  value,  Cd  —  1.338 /y/Re  =  0.013280.  Velocities  normalized  to  the  boundary  layer  edge  velocity,  uei 
versus  the  scaled  normal  coordinate,  7?,  are  shown  in  Fig.  9.4  for  the  two  finest  grids  in  one  family;  either 
computation  is  indistinguishable  from  the  Falkner-Skan  boundary  layer  analytic  result  that  accounts  for 
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FlG.  9.4.  Scaled  velocity  profiles  for  the  finite  plate  at  x  =  1.5;  Re  —  10,000;  g  —  (y/x) y/Rex /2. 


streamwise  pressure  gradient. 

The  largest  discretization  errors  as  well  as  the  largest  residuals  occur  near  the  leading  edge  singularity, 
as  can  be  noted  in  Fig.  9.2.  Although  not  tried,  a  local  refinement  near  this  boundary  would  be  beneficial. 

10.  Concluding  Remarks.  A  multigrid  method  for  solving  the  incompressible  Navier-Stokes  equa¬ 
tions  has  been  applied  to  a  classical  model  problem  of  fluid  dynamics:  flow  past  a  finite  fiat  plate  at  high 
Reynolds  number.  Elements  of  the  Full  Approximation  Scheme  multigrid  algorithm,  including  distributed 
relaxation,  defect  correction,  and  boundary  treatment,  have  been  presented  in  some  detail  for  the  three  main 
physical  aspects  encountered  in  the  simulation:  entering  flow,  wake  flow,  and  boundary  layer  flow.  Textbook 
efficiencies,  i.e.,  reduction  of  algebraic  errors  below  discretization  errors  in  one  multigrid  cycle,  and  residual 
reduction  rates  approaching  the  value  expected  for  elliptic  equations  of  nearly  one  order  of  magnitude  per 
cycle,  are  attained  for  second  order  accurate  simulations  at  a  laminar  Reynolds  number  of  10,000. 
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Appendix  A.  Linearized  Euler  Equations. 

The  linearized  Euler  equations, 


Loo(q)  =  0,  (A.l) 

with  periodicity  in  the  ^-direction  over  a  finite  domain,  0  <  x  <  L,  are  considered,  where  q  represents  a 
perturbation  from  freestream  values.  The  convection  operator  is  assumed  to  be  constant  as 

Qo  =  +  tdy ,  (A. 2) 

where  t  =  Uoo/^oo  represents  the  incidence  of  the  freestream  flow  with  the  a; -axis.  The  boundary  conditions 
are  taken  as  prescribed  velocity  components  at  inflow  and  pressure  at  outflow, 

(  W  )x=0  =  (  Uo  )eiu,y, 

V  Vo 

(  P  )x~L  =  (  Pl  )eiu,y-  (A.3) 

Brandt  and  Yavneh[5]  considered  entering  flow  (L  — >  oo)  with  inclusion  of  the  first  differential  approximations 
of  the  discrete  equations  to  confirm  algebraic  convergence  below  discretization  error  in  one  FMG  cycle, 
neglecting  boundary  effects.  Here,  we  consider  only  the  differential  solution  using  distributed  relaxation, 
q  =  MqqW,  and  include  boundary  effects.  Considering  w  of  the  form 


/  <*  \ 

|  b  I  e~axeiuy,  (A .4) 


then  Lqo  Mqo w  =  0  implies 


T  —a  +  icot 

0 

0  1  /a\ 

1  ° 

—a  +  uot 

0  1  1  6  1=0. 

(A.5) 

E  -a 

iu 

-a2  +  w2  1  '  c  ' 
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A  non-trivial  solution  (zero  determinant)  exists  for  values  of  a  as  below, 


a  C  iut)  a;,  —a;}. 


Thus,  the  general  solution,  w  =  wetoy,  can  be  written  as 


/  0  \  (  1  \ 

I  1  I  „-ivtx  d  ,  n  I 


1  1  1 

1  e~iutx  +  B2  I  0  1  e~iutx 

l  -i  1 

l  it  1 

w(l+t2) 

w(l+«2) 

/  0  \ 

/  0  \ 

I  0  I  e~ux 

+  B4  1  0  1  eux, 

1  1  ] 

1  1  1 

(A-6) 


(A.7) 


which  requires  four  boundary  conditions  to  close  the  system,  instead  of  the  three  required  with  the  prim¬ 
itive  equations.  Applying  a  Dirichlet  condition  for  at  inflow  supplemented  with  the  original  boundary 
conditions,  as  below, 


t  wi  ~  dx(w3)  \ 

/  «0  \ 

1  W2  1  = 

I  Vo  I 

s 

CO 

li 

o 

'  0  1 

(  Q0W3  )x=L  = 

(PL)eiuy, 

(A.8) 

the  coefficients  B\  —  B4  can  be  determined  and  are  given  below: 


Bi  = 

B2  =  uq  —  f(£uo  —  vq)ID2  +  2 plg  , 

Vi 

B3  =  ^TT)l{tu°-Vo)-ipLe~uL]’ 

Bi  =  M*«o  -  v0)e~2uL  -  ipLe~wL ] ,  (A.9) 

where  D\  =  l  +  e~2uL  and  £>2  =  1-  e~2uL.  Note  that  w2  is  a  function  of  vq  only;  w\  is  primarily  a  function 
of  uo  but  is  coupled  to  vq  and  pl  through  the  boundary  conditions,  Eq.  (A. 8);  the  coupling  is  rather  weak, 
however,  as  it  disappears  completely  for  vq  =  tuo,  as  is  usually  the  case,  and  L  — >  oo. 

The  primitive  variables,  q  =  qelujy ,  can  be  determined  from  q  =  M^w,  as  below, 


(A. 10) 
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(A.ll) 


Ai  =  uQ  +  tv o  +  y~-[i(tu0  -  vo)  +  2pLe  wL], 

U\ 

A2  =  +z(fu0  -  ^o)  +PLe~wL, 

As  =  -i(£u0  ^  v0)e~2uL  +Pl£~~u’Li 

It  can  be  verified  that  the  solution  above  satisfies  (A.l)  and  the  boundary  conditions  (A.3).  The  boundary 
conditions  for  w  discussed  in  the  main  body  of  the  text  are  discrete  forms  of  the  differential  boundary 
conditions  given  by  (A. 8)  above.  For  the  linearized,  constant  coefficient  case  considered  here,  both  the 
discrete  and  differential  forms  share  the  property  that  a  solution  to  the  distributed  relaxation  equations 
with  boundary  conditions  (A.8)  satisfy  identically  the  differential  equations  (A.l)  with  boundary  conditions 
(A.3). 

Appendix  B.  Tange ncy. 

The  linearized  Euler  equations,  Eq.  (A.l),  are  again  considered  but  with  t  =  0  (Qo  =  dx )  and  with 
periodicity  in  the  rr- direction  over  a  finite  domain,  0  <  y  <  H.  Linearized  tangency  boundary  conditions 
are  prescribed  as  v  at  the  top  and  bottom  of  the  channel, 


Considering  w  of  the  form 


(  v(x,0)  \  _(  v0  \ 
v(x,H)  vh 


f  a  \ 

|  b  I  e~ayeiux, 


then  Loo  Moo  w  —  0  implies 


T  icv 

0 

0  1 

f  a\ 

1  0 

iu 

0  1 

1  6  1 

iu) 

-a 

-a2  +  u2 

'  c  1 

A  non-trivial  solution  (zero  determinant)  exists  for  values  of  a  C  {a;,  —u}. 
Thus,  the  general  solution,  w  =  weiu}X,  can  be  written  as 


(B.l) 


(B-2) 


(B.3) 


/  0  \  /  0  \ 

w  =  Bi  |  0  j  euy  +  B2  |  0  j  e-uy,  (B.4) 

which  only  requires  two  boundary  conditions,  consistent  with  the  primitive  equations;  no  boundary  conditions 
can  be  given  for  wi  and  w2.  Applying  a  Neumann  condition  for  u>3,  as  below, 


(  (-9j,(w3))(x,0)  \  _(  v0  \  eit 
(-dy(w3)){x,H)  VH 
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the  coefficients  Bi  and  B2  can  be  determined  and  are  given  below: 

Bl  =  uD  [VH  ~  W°e_WHJ  ’ 

B2  =  cuD  _  VHe~UH ]  >  (B-6) 

where  D  =  1  +  e~2ujH . 

The  primitive  variables,  q  =  q_eiu)X,  can  be  determined  from  q  =  M^w,  as  below, 


/  -iu  \ 

(  -iu  ^ 

q  =  1  -w 

ewy  +  Bi 

- 

\  ioj  ) 

\  / 

(B.T) 


It  can  be  verified  that  the  solution  above  satisfies  (A.l)  and  the  boundary  conditions  (B.l),  recovering  the 
classical  aerodynamic  model  problem  for  the  flow  past  a  wavy  wall.  It  is  clear  that  ws  takes  the  role  of  the 
perturbation  potential;  the  Neumann  boundary  conditions  for  W3  are  implemented  discretely  at  the  tangency 
surfaces  in  the  main  body  of  the  text. 


Appendix  C.  No-Slip  Boundary. 

The  no-slip  boundary  conditions  along  the  plate  in  terms  of  the  ghost  variables  are 


[8w2](xc,  0)  =  [a'l(5w)3)](xc,0),  (C.l) 

=  [^(<5W3)](%,0),  (C.2) 

where  xc  denotes  the  x  position  of  the  cell-center  for  cell  (j.  k )  and  xg  =  xc  +  as  in  Fig.  C.l.  Since  a 
third  boundary  condition  for  the  ghost  variables  is  required  at  the  plate,  we  choose  to  split  Eq.  (C.l)  into 
two  separate  equations  as 
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[<5w2](a;c,0)  =  0,  (C.3) 

[a?(««*8)](*e,0)-0.  (C.4) 

At  the  location  (zfl,0),  Swi  and  SW3  can  be  approximated  in  terms  of  nearby  values  as 

[5u»i](a:fl,0)  =  ^(<5wi)j,2  +  ^(<5wi)j,i,  (C.5) 

=^[(^3)j+i,2  -  (Sw3)jt2  (C.6) 

+  (<5w3)i+i,i  -  (<5w3)j,i], 


In  relaxing  the  jth  column,  we  assume  that  =  0.  Prom  Eq.  (C.4),  we  also  have  (Sws)jti  —  (<ftu3)j,2- 

Then  Eq.  (C.2)  can  be  written  as 

=  -(^l)j,2  -  ^(<^3)^2,  (C.7) 

which  is  an  implicit  boundary  condition  equation  to  be  implemented  in  relaxing  Eq.  (4.6)  at  the  wall. 

Now  assume  the  convection-diffusion  operator  is  constant,  defined  with  a  computational  stencil  as  below, 

r  cN  1 

Qv  =  I  cw  co  ce  I  •  (C.8) 

cs 

In  a  lexicographic  pointwise  relaxation,  the  matrix  to  solve  for  the  (Jw)^  values  is  as  below, 

r  Co  0  0  1/  Sw2  \  /  r2  \ 

I  0  co  —  cs  —  2csh~1  II  6w\  I  =-|  ri  I  •  (C.9) 

h~l  h -1  3  h~2  Sw3  T3 

3, 2  3, 2 

This  system  couples  the  implicit  equations  for  Swi  and  Sw^  at  the  cell  adjacent  to  the  no-slip  boundary, 
necessitating  a  local  2x2  block  matrix  solution.  After  solving  for  (£w)J}2  (and  thereby  (Sws)jti  and  (6iui)j,i) 
and  changing  the  primitive  variables  through  M£w,  it  can  be  shown  that  the  updated  residuals  of  cell  (j,  k) 
are  zero.  For  variable  coefficients  in  the  convection-diffusion  operator,  the  residuals  differ  from  zero,  as  they 
do  in  the  interior  of  the  mesh. 

Considering  relaxation  of  the  entire  column  of  cells,  the  implicit  equations  for  the  cells  away  from  the 
boundary  remain  in  lower  triangular  form.  Thus,  the  equations  can  be  solved  using  an  LU  decomposition 
with  only  a  small  overhead.  In  this  instance,  the  entire  column  of  residuals  are  zeroed  out  for  a  constant 
coefficient  convection-diffusion  operator. 

Appendix  D.  Boundary  Stencil  Modifications. 

The  four-point  upwind-biased  stencil  considered  here  requires  special  treatment  near  boundaries.  For 
prescribed  velocity  boundary  conditions  at  inflow,  a  modification  is  required  at  j  =  2  for  the  x— momentum 
equation  and  at  j  =  2  and  j  =  3  for  the  y -momentum  equation.  For  prescribed  pressure  at  outflow,  a 
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modification  is  required  at  j  —  Nx  for  both  momentum  equations.  The  simplest  approach  to  maintain 
second  order  accuracy,  used  herein  for  the  wake  and  boundary  layer  simulations,  is  to  use  a  first-order  two- 
point  stencil  for  these  points  near  inflow  and  a  fully-upwind  stencil  (k  =  -1)  near  outflow.  For  the  entering 
flow  simulation,  more  accurate  stencils  were  used  at  inflow  as  shown  below;  u^k  and  vz/2,k  represent  given 
boundary  values  at  x  —  0,  as  in  Fig.  6.1. 

Considering  the  x-momentum  equation,  for  the  SUD-2  scheme,  the  required  term  ux  is  computed  using 
nearby  points  and  the  gradient  at  x  =  0,  i.e., 


(dxu)\h,y  =^[-5«l,fc  +  4ti2  ,fc  +  U3,fc] 

-\{dhxu)\o,y,  (D.l) 

where  (h,y)  denotes  the  vertical  interface  midpoint  of  the  (2 ,fc)  cell  and  (d^u) |o,y  =  ~(dyv) |o,y  is  given  at 
inflow  from  continuity,  as 


(flyu)kv  =  24 ^[27v3/2,fc  27v3/2jfc-i 

“  ^3/2, fc+1  +  ^3/2,/c— 2]  •  (D.2) 

For  the  NUD-2  scheme,  central  differencing  (k  =  +1)  is  used. 

Considering  the  y— momentum  equation  with  either  scheme,  central  differencing  is  used  at  the  j  —  2 
column  of  cells  and  a  third-order  4-point  formula  at  the  j  =  3  column  of  cells,  i.e., 

(9xv)\h/2  ,y  =  ^[-4vs/2,fc  +  3^2,  A:  +  (D‘3) 


(dxV)Uh/2,y  =  ^[16«3/2,fc  -  45 V2,k 

+  20vZ)k  4*  9v4 ,/s],  (D*4) 

where  (h/2,y)  and  (3h/2,y)  denote  the  horizontal  interface  midpoints  of  the  (2,  &)  and  (3,fc)  cells,  respec¬ 
tively. 
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